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Summary: The Lagrangian perturbation theory on Priedman-Lemaitre 
cosmologies investigated and solved up to the second order in earlier pa- 
pers (Buchert 1992, Buchert & Ehlers 1993) is evaluated up to the third 
order. On its basis a model for non-linear clustering applicable to the 
modeling of large-scale structure in the Universe for generic initial condi- 
tions is formulated. A truncated model is proposed which represents the 
"main body" of the perturbation sequence in the early non-linear regime 
by neglecting all gravitational sources which describe interaction of the 
perturbations. However, I also give the irrotational solutions generated 
by the interaction terms to the third order, which induce vorticity in 
Lagrangian space. The consequences and applicability of the solutions 
are put into perspective. In particular, the model presented enables the 
study of previrialization effects in gravitational clustering and the onset of 
non-dissipative gravitational turbulence within the cluster environment. 



1 



1. Introduction 



For a long period of research in cosmology inhomogeneities in the Universe have 
been modeled on the basis of a perturbative approach exploiting the instability 
of the standard cosmologies of Friedman-Lemaitre type against perturbations of 
the density and the velocity field. This approach is Eulerian, i.e., the perturba- 
tions are evaluated as a function of Eulerian coordinates (see, e.g., Peebles 1980). 
The limitations of this approach have been widely recognized, since it relies on 
the smallness of physical densities which is inappropriate for the modeling of the 
high density excesses observed in the Universe. Zel'dovich (1970, 1973) has real- 
ized this situation and has proposed an approximate extrapolation of the linear 
perturbation solution into the non-linear regime. 

In recent papers the Lagrangian theory of gravitational instability of cosmolo- 
gies of Friedman-Lemaitre type for preasureless ("dust") matter has been in- 
vestigated and solved up to the second order in the deviations from homogene- 
ity (Buchert 1992, Buchert & Ehlers 1993, henceforth abbreviated by B92 and 
BE93). This theory does not rely on the smallness of the density of the inhomo- 
geneities, only the deviations of the particle trajectories from the homogeneous 
Hubble flow are treated perturbatively. This is possible, since the field of trajec- 
tories f{X, t) is the only dynamical variable in the Lagrangian picture. Interest- 
ingly, the widely applied "Zel'dovich approximation" for modeling the formation 
of large-scale structure in the Universe was found to be contained in a subclass 
of first order irrotational perturbation solutions in this theory (B92). The first 
order solutions have been analyzed in an earlier paper (Buchert 1989) in which 
they are shown to provide exact three-dimensional solutions in the case of locally 
one-dimensional motion, i.e., when two eigenvalues of the peculiar- velocity gradi- 
ent vanish along trajectories of fluid elements. The general plane- symmetric case 
{globally one-dimensional motion) constitutes a subclass of these solutions. 

It should be emphasized that solutions of the Euler-Poisson system, exact 
or in a perturbative sense, will depend non-locally on the initial conditions for 
the inhomogeneities. In general, they have to be constructed by solving ellip- 
tic boundary value problems, as will be discussed in full detail. Despite this, 
Zel'dovich's approximation is local in this sense. It can be made rigorous that 
the first order perturbation solutions can be made local without loss of generality 
(B92). This will be recalled in APPENDIX B. The problem of uniqueness in 
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general perturbation solutions will be considered in a separate paper (Ehlers & 
Buchert 1993). 

The success of Zel'dovich's model as an approximation for irrotational self- 
gravitating "dust" flows in the early non-linear regime is commonly appreciated. 
Its range of applicability can be roughly limited to the epoch shortly after the 
first shell-crossing singularities in the flow develop, and provides an excellent 
approximation for the density field down to the non-linearity scale (i.e., where 
the r.m.s. deviations of the density from homogeneity exceed unity) in com- 
parison with numerical N-body simulations (Coles et al. 1993). In contrast to 
other analytical models for the formation of large-scale structure, the Lagrangian 
perturbation theory offers a systematic and explicit way to extend the range of 
validity of Zel'dovich's model for generic initial conditions (Buchert 1993). While 
the first order approximation chiefly covers the (up to the epoch of shell-crossing 
dominating) kinematical aspects of the structure formation process, the second 
order approximation firstly involves the tidal action of the gravitational field. The 
collapse process is significantly accelerated by this action. Also, tidal forces con- 
stitute the essence of so-called previrialization effects in gravitational clustering 
(Peebles 1990, BE93). 

Zel'dovich's model has shortcomings after the epoch when shell-crossing sin- 
gularities develop, since shells of matter are predicted to cross freely, and the 
kinetic energy of particles at the caustic exceeds the gravitational potential en- 
ergy. As a result no potential well is formed in contradiction with numerical 
simulations. The advantage of high-spatial resolution N-body simulations is obvi- 
ous especially with respect to the modeling of multi-stream systems which develop 
after shell-crossing. Analytically, this modeling is not straightforward, since the 
particles are not treated as interacting bodies. Rather, the particles are viewed 
as tracers of a flow, which generally looses uniqueness as shells of matter over- 
lap. However, while the Eulerian representation of any solution breaks down as 
soon as singularities in the density field (caustics) form, the Lagrangian represen- 
tation allows to follow the trajectories / of fluid elements across caustics. The 
solution / remains regular, only the transformation from Eulerian to Lagrangian 
space is singular. Non-uniqueness of the flow is only realized in Eulerian space, 
where a patch of matter can originate from different Lagrangian 'particles', we 
say that the flow consists of several streams in such a region. Therefore, analyti- 
cal Lagrangian approximations are in principle capable to model fully developed 
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non-linear situations including multiple shell-crossings. 

The phenomenology of a self-gravitating multi-stream system is complex, the 
number of streams is systematically increasing in time. Numerical simulations 
demonstrate that a hierarchy of singularities forms (e.g. Doroshkevich et al. 1980), 
shells do not cross freely without limit as e.g. predicted by Zel'dovich's model, 
rather the gravitational action of the firstly formed three-stream system (called 
'pancakes' in the cosmological context) forces the inner fluid elements to recollapse 
and participate in a second shock structure and so on. This way a five-, seven-, 
nine-stream system etc. develops. A potential well is established leading to a self- 
trapped quasi-stationary structure around the firstly developed singularities. The 
morphogenesis of such a non-dissipative structure (e.g. a cluster of galaxies in the 
present context) is similar to that of a flower: the singular density peaks split up 
into shock structures which subsequently emerge from and move away from the 
center. Apparently, a cluster is composed of expanding shells, while permanently 
new shells are created in the center. A nested moving system of singularities is 
formed on smaller and smaller spatial scales. In real terms the situation is more 
complicated: in addition higher-order singularities (such as the 'breaking' of a 
pancake boundary) increase the number of streams as was originally discussed 
for Zel'dovich-Arnol'd pancakes by virtue of analyzing a numerical simulation in 
two spatial dimensions (Arnol'd et al. 1982). The resolution of higher-order sin- 
gularities can also be appreciated in high-resolution plots of analytical mappings 
such as the "Zel'dovich approximation" (in 2D: Buchert 1989, in 3D: Buchert 
Bartelmann 1991). 

Analytical models have also been used to analyze self-gravitating multi-stream 
systems: Approximating a generic density peak locally by a triaxial ellipsoid, a 
model for the evolution of this peak has been formulated by Gurevich and Zybin 
(1988a,b) by performing a transformation of the known spherically symmetric 
solution. They show that the density distribution around the collapsed peak can 
be represented by a scaling law with a power of —1.8 in this model. This law 
is found to hold on a wide range of spatial scales from the scale of clusters down 
to the scale of typical star distances. Also Moutarde et al. (1991) found a similar 
scaling behavior of the density field in the early non-linear regime of the collapse. 
The correspondence to the scaling behavior of observed galaxy concentrations 
as measured by the two-point correlation function is striking and suggests that 
multi-stream hierarchies in the density field of dark matter might be responsible 
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for the clustering properties of galaxies (Gurevich & Zybin 1988a, Berezinsky et 
al. 1992). However, this conclusion depends on the assumption that the dark 
matter dominates the matter density down to these scales (Gurevich & Zybin 
1990). It is likely that the fraction of baryonic matter depends on scale and posi- 
tion which implies that the influence of hydrodynamic or radiative effects of the 
baryonic component might alter this picture. Also, it is not shown that a generic 
collapse occurs at the peaks of the inzim/ density fleld. It is likely that, in general, 
a collapse occurs closer to the trajectory of the maximum of the largest eigenvalue 
function of the initial displacement tensor (Shandarin, priv. comm.). Gurevich 
and Zybin call this hierarchical structure "non-dissipative turbulence". Its exis- 
tence is related to the early conjecture by Mandelbrot (1976) that singularities in 
self-gravitating flows might form a fractal, although they rather display a multi- 
fractal scaling (Martmez 1991, Jones 1993). Indeed, a self-similarity of caustic 
patterns can be appreciated (see BE93) in accord with the existence of self-similar 
solutions (Filmore h Goldreich 1984, Bertschinger 1985). Moreover, in developed 
gravitational turbulence the density scaling law seems to be generic, i.e., it seems 
not to depend on the initial fluctuation spectrum according to Gurevich and Zy- 
bin's model. However, it appears, at least with the amount of non-linearity one 
can resolve, that the impact of large-scale power will destroy the convergence 
to this local power law form, thus retaining the memory from the initial condi- 
tions, (compare the flnal power spectra for difl"erent intial conditions in Melott &; 
Shandarin 1990). 

The value of Lagrangian perturbation solutions has to be tested against the 
phenomenology of "non-dissipative turbulence" . Indeed it has been demonstrated 
that the second order solutions do describe the onset of such a hierarchy. They 
predict a second shell-crossing singularity within Zel'dovich-Arnol'd pancakes. A 
second bifurcation branch appears on the critical manifold of the flow in La- 
grangian space (BE93). This prediction suggests that m-th order Lagrangian 
perturbations will transform pancakes into 2m-|-l-stream systems in the coarse 
of time. The stage until when the perturbation solutions are valid could be esti- 
mated by the time when the shell-crossing singularities of the corresponding order 
appear. 

In the present paper, I evaluate the third order perturbation solutions. In 
this line the need for higher order solutions for the purpose of modeling large- 
scale structure should be questioned. Certainly, it will not be adequate to derive 
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higher orders unless we can expect to obtain an appropriate model for large-scale 
structure which remains sufficiently simple to handle with respect to applications, 
apart from the fact that the derivation of higher order perturbation solutions for 
generic initial conditions is cumbersome. As an argument to go to the third 
order I consider the structure of the equations to be solved: they are cubic in 
the basic dynamical variable (see B92, APPENDIX), so we can expect that a 
third order solution will cover the main effects of the perturbation sequence in 
the early non-linear regime. We first derive the longitudinal perturbations, i.e., 
for the case where the perturbations admit a potential in Lagrangian space. In 
contrast to opinions stated in the literature (Moutarde et al. 1991, Bouchet et al. 
1992, Lachieze-Rey 1993), this restricts the generality of the solutions as will be 
discussed in detail. The restriction vanishes if interaction of perturbations (among 
the first and second order perturbations here) is neglected, an assumption which 
will define a class of third order models which we consider the "main body" of the 
perturbation sequence. Since, in general, the interaction terms generate vorticity 
in Lagrangian space (even for irrotational ^ows considered throughout this paper), 
there will be a transverse part of the third order solution, which I shall derive. 

As an example the approximation is evaluated on an Einstein-de Sitter back- 
ground for initial conditions which correspond to Zel'dovich's model, and the 
result is expressed in terms of the initial conditions for the peculiar-velocity po- 
tential, or the peculiar-gravitational potential, respectively. With this specifica- 
tion of the initial conditions (discussed in section 3) we have assumed a functional 
relationship between the peculiar-velocity potential and the peculiar-gravitational 
potential (here: proportionality). Relaxing this relationship will also induce vor- 
ticity in Lagrangian space already at the second order level (see BE93). 
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2. The perturbation equations in Lagrangian form 



Let us briefly summarize the formal consequences of a Lagrangian description of 
self-gravitating flows in Newton's theory. We introduce integral curves x = f{X,t) 
of the velocity fleld v{x,t): 

%=v{f,t) , f{X,to)=:X . (1) 

These curves are labelled by the Lagrangian coordinates Xi ; Xi are non-rotating 
Eulerian coordinates. We can express all Eulerian fields in terms of the field of 
trajectories / as was explained in Buchert (1989), B92 and BE93. We denote the 
determinant of the deformation tensor (fi^k) by J, where the comma denotes par- 
tial differentiation with respect to the Lagrangian coordinates, a dot will denote 
Lagrangian time derivative ^ := dt\x + v - Vx — dt\x', comma and dot commute. 
Recall that mass conservation is guaranteed in the Lagrangian representation irre- 
spective of any equations the trajectories / might obey. The Euler-Poisson system 
for "dust" matter can be cast into a set of four evolution equations for the single 
dynamical variable / (Buchert Sz Gotz 1987, Buchert 1989): 

^Pib Q(^Xi X2 Xs) ^ ' (2a, 6, c) 

E V-^^ mx'x^X^ = -^^Gp{X) ; p(X)>0 . (2d,e) 

aTi 2 ^[X^,X2,X^) 

(indices run from 1 to 3, if not otherwise explicitly stated; henceforth, Vq denotes 
the nabla operator with respect to the Lagrangian frame which commutes with 
the dot). 

We proceed to the evaluation of the perturbation equations. (The reader may 
consult earlier papers for more details.) We first make the following perturba- 
tion ansatz for longitudinal perturbations superposed on an isotropic homogeneous 
deformation: 

/ = a{t)X + p ; p = eVoV'^'^ + VqV'^'^ + VqV'^^^ • (3a) 

The parameter e is supposed to be small and dimensionless. It can be considered 
as the amplitude of the initial perturbation field. We formally split the initial 
density accordingly: 

p = PH + sSp^'^ + sHp^^^ + eHp^^^ . (36) 
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However, we can choose initial data at our convenience, so we put 

5^(1) := Sp =: phS , := , ^p^^^ := , (3c) 

o 

where 5p denotes the total initial density perturbation, 5 the initial density con- 
trast. This formality is adequate, since the density needs not be perturbed in the 
Lagrangian framework. 

In what follows, /(V',i,fc) = tr('0,i,fe) = AqV^, II{^,i,k) = i[(^^'(V',z,fc))^ - ^^((-0,^,^)^)] 
and III{ip^i^k) = det('0,i,fe) denote the three principal scalar invariants of the tensor 
{ip,i,k)- For the derivation of the perturbation solutions we shall first concentrate 
on the source equation (2d) and consider the conditions (2a,b,c) as constraint 
equations, which are checked after the solution is obtained; the constraint equa- 
tions and the resulting constraints on initial conditions for the longitudinal third 
order approximation are given in APPENDIX A. We shall then construct the 
transverse part which removes the constraints obtained for the pure longitudinal 
solution. 

Inserting the ansatz (3) into equation (2d), we obtain the following set of equations 
to be solved: 



£°|-47rG'pij = Sda^ - o^a} = 



£i{-47rG5p(i) = [(2aa - a^A) + ^] AqV'^'^} = 
e^[-4nGSp^^^ = [{2da - a^A) + a" ^] AoV'^'^ 



(4a) 
(46) 



r- ,(l)^ V- 9(V' i\ V V , ^c) ^ 

a,b,c 

l L dt'^i 







(4c) 



+{a - aA) Y: ea,c q^x,,X,,Xs) + ^ ''''^ 

a,b,c 



a,b,c 



a(Xi,X2,X3) ^ a(Xi,X2,X3) 



+ 



X 9(x,,x2,x3) - ^ = 

a,b,c 



Each order (starting from e^) contains a term where a linear operator acts on the 
different potentials Similarily, a quadratic and a cubic operator acts 
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on quadratic or cubic invariants formed by second derivatives of However, 
at the third order, there appears a quadratic term describing the interaction of 
hnear and quadratic perturbations. Higher orders will subsequently add such 
interaction terms between the potentials at different or equal orders. 

We now derive the solutions of the equations (4) in the case of a "flat" back- 
ground. I present a simplifled derivation for a special class of initial conditions 
which is relevant for the modeling of large-scale structure. Therefore, it is illus- 
trative to derive all orders with the same restriction and procedure. The reader 
who is interested in applications of the solutions only may proceed directly to 
section 4. 
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3. Solution of the perturbation equations 
on an Einstein-de Sitter background 



In what follows we make use of a restriction of the initial conditions to simplify 
the derivations. We require that, initially, the peculiar-velocity u{X,to) be pro- 
portional to the peculiar-acceleration w{X,to): 

u{X,to) = w{X,to)to , (5) 

where we have defined the fields as usual (compare Peebles 1980, B92). This 
restriction has proved to be appropriate for the purpose of modeling large-scale 
structure since, for irrotational flows, the peculiar-velocity field tends to be par- 
allel to the gravitational peculiar-field strength after some time. The reason for 
this tendency is related to the existence of growing and decaying perturbations, 
the growing part supports the tendency to parallelity. The assumption of irrota- 
tionality should be adequate down to the non-linearity scale. However, besides 
the possibility of non-linearly enhanced primordial vorticity (B92), shell-crossings 
generate vorticity on scales below the non-linearity scale (see, e.g., Doroshkevich 
1973, Chernin 1993 and ref. therein). We should keep this in mind, since a third 
order approximation actually should give a proper description of these smaller 
scales. The treatment of vorticity is no longer academic, but might play a de- 
cisive role for the dynamics within clusters of galaxies. Also, decaying solutions 
should be considered with some care in the non-linear regime. They couple to 
the growing solution, the role of decaying and growing solutions in an expanding 
environment can be interchanged in a collapsing environment as was discussed by 
Gurevich and Zybin (1988b). We note that the restriction (5) is commonly used 
in the literature. The tendency of the flow expressed by (5) has been proved for 
general flrst order (Bildhauer &: Buchert 1991) and a large class of second order 
irrotational flows (BE93). Note that in the case (5) we have to give one initial po- 
tential only, whereas the general initial value setting for irrotational flows would 
require two. Any restriction of this type simplifies the calculations enormously, 
but it implies the restriction to irrotational fiows. (An alternative restriction of 
this type has been discussed in B92 and BE93). 

We shall use in the following the initial peculiar-velocity potential S defined 
as u{X,to) =: VoS{X). The initial peculiar-gravitational potential 0, w{X,to) =: 
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— Vo0(-^) is related to it as <S = 
of the equations (4) of the form: 



(j)to (eq. (5)). Consequently, we seek solutions 



(6a) 
(66) 
(6c) 



where the potentials 5^^^ and S^^^ have to be related to the initial condition S. 
We shall see that this relation will require the solution of elliptic boundary value 
problems expressing non- locality of the solutions. Formally we require qzito) — 0, 

qzz{to) = 0, Qzzzito) = 0. 

3.1. The zero order solution 

One class of solutions of the Euler-Poisson system (2) is formed by the homoge- 
neous and isotropic Friedman-Lemaitre cosmologies: 



equation (4a). Its general solution is given by solutions of Friedman's differential 
equation as an integral of (4a): 

+ const SttGoh + A 
= 3 ' 

where pn = PHa~^ is the background density. Henceforth we restrict all consid- 
erations to the Einstein-de Sitter case {const = 0,A = O). Then, the zero order 
solution reads: 



For convenience, we shall express all time-dependent coefficients in terms of the 
solution a. The constant —A^nGpu will be written in terms of its value —i^ for 



the solution (7c). 

3.2. The first order solution 

By virtue of the ansatz (6a), equation (4b) simplifies to the following equation 
(in the sequel we put A = and use (3c) and (7c)): 





(7c) 




(8) 
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Using Poisson's equation for the initial potential and the relation S = —(f)to we 
have: 

^oS=-^S . (8a) 

Hence, solutions of (8) can be found as solutions of the linear ordinary differential 
equation: 



a 1 
with: 

Ao<S(^) = Ao5 to = i,fc) ^0 • (8c) 

The solution to (8b) consists of two linearly independent solutions of the homo- 
geneous part: 

,W = C^a^ , = C(^)a- , (9a) 

and a particular solution: 



(96) 



The coefficient C^^^ is found by inserting q^^ into (8b); the coefficients C^^ and 
6*2^^ are found by the requirement that the coefficient functions of the peculiar- 
velocity and -acceleration equal to 1 at t = to for the solution q^^^ = g^^^ + ?2^^ + • 
Restricting the general solution q'^-^^ according to (5) we find (B92): 

q, = ^{a^-a) . (10) 

A discussion of the uniqueness of this solution related to the use of Sto instead of 
5(1) is to be found in APPENDIX B. 



3.3. The second order solution 

Inserting the ansatz (6b) into the equation (4c) we find: 

(,„ + 2%„)A„.<^. = -(M.4iL)2n(.S). (11) 

Solutions to (11) can be found by solving the linear ordinary differential equation 
(inserting the first order solution g^): 

q^^ + = - ^ ——\-.Q^ \t) , (11a) 

a \ a I a'^ J 
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with: 

Ao5(^) = 2II{S';11) ■ (11^) 

The solution to (11a) consists of two hnearly independent solutions of the homo- 
geneous part: 

,f ) = , = )a-i , (12a) 

and a particular solution: 

i'^ = C^'^ (gf^ I'gf ^(^)dt-gf j\[^^gi^}d?j=C(^\-l + la-') . {12b) 

The coefficient C^^^ is found by inserting q^^ into (11a); the coefficients C^^^ and 
C2^^ are found by the requirement that the coefficient functions of the pecuHar- 
velocity and -acceleration vanish at t = to for the solution q^^^ = q^"^ + q^^ + q^'' . 
We find for the restriction (5) the following result (compare BE93, section 5): 

/3\% 3 o 3 2 1 4 1, 
^-= 2 ^-14^+5^-2^+35^"') • 



3.4. The third order solution - longitudinal part 
Inserting the ansatz (6c) into the equation (4d) we find: 

- (^^'^^O • (14) 

To obtain solutions of (14) we use the linearity of Poisson's equation and spHt the 
potential S^^^ into a part S^^"'^ generating the cubic source term, and a part 5^^'') 
generating the quadratic source term of the interaction of first and second order 
perturbations. We then have to solve separately the Unear ordinary differential 
equation (inserting the first order solution g^): 

Qzzz + 2-g... = -^(izQl =: S^^^Ht) , (14a) 

ft ft 

with: 

AoS^''^^ = 3III{S%) , (146) 
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and the linear ordinary differential equation (inserting the first and second order 
solutions Qz and Qzz)'- 



qzzz + '^-qzzz = --^qzzqz + -{qzzqz + qzqzz)=-Q^^''\t) , (I4c) 

a a'' a 

with: 

A general solution to (14a) or (14c), respectively, consists of two linearly inde- 
pendent solutions of the homogeneous part: 

^(3a,6) ^ ^(3a,6)^2 ^ ^(3a,6) ^ ^(3a,6)^-i ^ ^^^^^^^ 

and a particular solution for each source term: 

qi^Sa)^C^Sa) |^^(3a) g(3«)^(3a)^^ _ ^(3a) gi^a) 

= C^3b)\^-a-- + -a-^--a-i) . (166) 
4^7 10 2 35 ^ ^ ^ 

The coefficients C*^^"'''^ are found by inserting q^"'"^^ into (16a) or (16b), respec- 
tively; the coefficients C^°"^^ and c'^"''^^ are found as in the second order case 
by the requirement that the coefficients of the initial peculiar-velocity and - 
acceleration vanish for the general solutions. I find for the restriction (5) the 
following result: 

„ /3\^, 1480 Sol 16 _i. . 

^ fsY , 5 4 33 3 7 2 1 4 1 4 _i, 

ql^^^[-\ ( — a* a-^ ^ a -\ -\ a 2) . (176) 

'izzz \^ 70 10 2 35 105 ^ ^ ' 
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3.5. The third order solution - transverse part 
Inserting the longitudinal ansatz 



|2 



= Vo ( + g,5(i) + g,,5(2) + g«,,5(3«) + q^^^S^^'^ \ (18) 



into the equations (2a,b,c), we find restrictions on the initial conditions. These 
restrictions are derived in APPENDIX A. Note that, instead of the equations 
(2a,b,c) for the irrotationality of the gravitational field strength, we can use the 
corresponding equations for the irrotationality of the velocity in the case of irro- 
tational flows (see the LEMMA proved in BE93), which we shall do. 
The restrictions only arise at the third order level. This implies that the third 
order solution is not purely longitudinal, it is not fully covered by the ansatz (18). 
In the following we seek a transverse part of the third order solution by extending 
the ansatz (18): 

f = f^ + F ; F--=qL.^ ; S:=-Vox5(3-) . (19a) 

We have introduced the vector potential S^^^\ on which we impose the following 
gauge condition (compare APPENDIX B): 

Vo-.S^^''^ . (196) 

With the ansatz (19a) we generate no additional equations to be fulfilled in the 
source equation (2d). The integrability conditions for the velocity (BE93: equa- 
tions (5d,e,f)) only yield an equation of the order as expected {i,j,k = 1,2,3 
cyclic ordering): 



e'[a {aql,, - aql,,) (Vo x (-Vq x S^^^^^)) ^ 

+a{qzzqz-qzqzz)epq[j-^^j^-^-—yj =0 . (20a, 6, c) 

Using the vector identity Vo x (Vo x S^^"^) = Vo(Vo • S^^''^) - AqS^^^^ and (19b), 
we find solutions of the equations (20) by solving the linear ordinary differential 
equation (inserting the first and second order solutions q^ and g^^): 

-q%, = -- {q,,q, - q,q,,) =: g^^'^\t) , (21a) 



izzz 

a a 
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with: 

:(2) c(i) 



A general solution of (21a) is given by: 



with the integrating factor M{t) = e -^^o " . The coefficient C^^^^ is found by the 
requirement q^^^{to) — 0. We finally obtain: 

/SV/ 1 4 3 3 1 2 1 4 1 4 _i. 
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4. Result and Discussion 



4.1. The solution 



THEOREM 



With a superposition ansatz for Lagrangian perturbations of an Einstein- de Sitter 
background we have obtained the following family of trajectories x = f{X,a) as 
irrotational solution of the Euler-Poisson system up to the third order in the per- 
turbations from homogeneity. The general set of initial conditions {(j){X),S{X)) is 
restricted according to <S = — 0to- (The parameter e is considered as the amphtude 
of the initial fluctuation fleld; a{t) = (t/to)'^/^, i,j,k = 1,2,3 with cyclic ordering): 

f = aX + q,{a)VoS^'\X) + q,,{a) VoS^^\X) 

+ ql^M WoS^'-\X) + Vo«S(^^)(X) - q^^^^ia) Vo X S^^'^\X) , (24) 



with: 



and: 



(a^-a) , (24a) 



33301 4 1. , 

2) (-14" +8" -2"+ 35°"') ' (2^^) 

n /3\'^,l/l3q39l 16 1, 

= ( 2 j (-9" + r - r + r - ir^"'''^ • ^^^^ 

u / 3\^ , 5 4 33 o 7 2 1 4 1 4 1, „ 

''-=^2) ^42^^ -70^ +10'^ - 35^^ + 105^"^^ ' ^'''^ 

/ 3\% 1 4 3 3 1 2 1 4 1 4 _i, , 

_ ( — a-^H + -a 02 H a 2) (24e) 

Hzzz I 2/ ^14 14 10 2 7 35 ^ ^ ^ 



Ao^W = 1(5 ,,fc) to , (24/) 

Ao5(2) = 2//(5« ) , (24^7) 

Ao5(3«) = 3///(4'l) , (24/i) 

a,o,c 
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REMARKS: 

The potential ^S^^**) and the vector potential .S^^'^^ generate interaction among the 
first and second order perturbations. The general interaction term is not purely 
longitudinal. In order to satisfy the Euler-Poisson system with the longitudinal 
part only, we have to respect the following constraint (APPENDIX A): 

Vo5^2) ^ yy^y^^W) _ (24m) 

The potential 5*^'^^^ generates the symmetric part, whereas the vector potential 
5(3c) generates the anti-symmetric part of the interaction. 

Proof: The proof follows by inserting the solution (24) into the perturbation equa- 
tions (4), which has been done using the algebraic manipulation system REDUCE. 

We now discuss some issues which are related to the practical use of the solution 
(24) as a model for the evolution of large-scale structure. 

4.2. The construction of local forms 

The usefulness of Zoca/ approximations has been pointed out by Nusser et al. (1991) 
who apply the "Zel'dovich approximation" as a tool for locally reconstructing 
the density field from observed peculiar-velocity data. Contrary, solutions of 
Newtonian equations for the evolution of self-gravitating "dust" continua are non- 
local., since, e.g., the gravitational potential has to be a solution of the (elliptic) 
Poisson equation. This solution involves integrals over large space regions. As 
far as the first order Lagrangian perturbation solution is concerned, this is not a 
contradiction: without loss of generality, i.e., without change of physical quantities 
like the density contrast or the divergence of the peculiar- velocity, we can use the 
initial condition Sto instead of (APPENDIX B). With this reduction the 
iterative procedure of constructing the displacement vectors in (24) simplifies: 
the source terms are (except the interaction terms (24i-l)) completely expressible 
in terms of second derivatives of the initial potential S. 

The situation is more delicate at higher order levels. We no longer are able 
to write second or higher order terms as a local approximation. Although it is 
computationally simple to obtain the displacement vectors in (24) by employing 
Poisson solvers (see 4.3), it is useful to know which classes of initial conditions 
admit the construction of local forms. With such closed form expressions the 
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study of special solutions is explicitly feasible without using a numerical Poisson 
solver. 

According to COROLLARY 1 proved in (BE93), a local form can be obtained 
for second order displacements. It reads: 

Vo5(2) = VqS (AqS) - {VqS ■ Vo) VqS ; Vo<S x Aq Vq^ = . (25a, b, c, d) 

The local form (25a) is constructed such that its divergence agrees with the source 
term in (24g), its curl is, however, in general non-zero, it only vanishes if (25b-d) 
is statisfied. The latter equations express the fact that the form (25a) cannot be 
used in general, constraints have to be obeyed. 

Similarily, we can ask for a local vector form whose divergence agrees with the 
third principal scalar invariant of (5i,fc), which is the source term in equation 
(24h). Indeed, the following expression has the required property: 

COROLLARY 1 

The vector Vo«S*^^") with the components 

{VoS('^\ = J2i^oS^'^),jf,, (26a) 

i 

has the property: 

Ao5(3«) = 3m(5W ) , 

where Jf/^ are the subdeterminants of the tensor {S^l\). The following constraints 
have to be satisfied in order that Vo^^"^"^ is curl-free: 

J](Vo5(i)),iJ5,,,.]=0 , k^j . (266, c,d) 

i 

(The proof is done by explicit verification). 

The source term in (24i) which describes the longitudinal part of the interaction 
of first and second order perturbations has a similar structure as the second order 
source term. We can construct a local form by analogy. The integral of the 
interaction source term in equation (24i) reads (BE93, COROLLARY 1): 

Ai (Vo5(2)(Ao5(^)) - (Vo<S(2) • Vo)Vo5^^))+A2 (v qS^^^AqS^''^ - Vo^^^^ • Vo)Vo5(')) 

(27a) 
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We have taken the hnear combination of the two possible integrals as a general 
integral, where we have to assure Ai + A2 = 1. In order to satisfy the requirement 
that the vector (27a) be a solution of the Poisson equation (24i), we have to assure 
that it is curl-free which implies (BE93, COROLLARY 1): 

Ai (VoS^^^ X AoVo5^i)) + A2 (Vo^^i) x AqVqS^^^^ = . {27b, c,d) 

We finally give an integral expression for the source terms in (24j,k,l) which de- 
scribe the transverse part of the interaction of first and second order perturbations: 

COROLLARY 2 

The vector S: 

S := -Vo X S^^""^ = 
f^i ((Vo5(2) . ^^^^^^(1)^ ^ (-(Vo^^i) • Vo) Vo5(2)^ (28a, b, c) 

has the property: 

(Vo X ^). = ,^^J ■ 

We have taken the linear combination of the two possible integrals as a general 
integral, where we have to assure fii + fi2 = ^- In order to satisfy the require- 
ment that the vector components (28a, b,c) be solutions of the Poisson equations 
(24j,k,l), we have to assure that the displacement vector S can be represented in 
terms of the vector potential S^^'^^: S = — Vq x S^^'^\ i.e., the vector field S has to 
be source-free: 

/^iVo • ((Vo5(2) . ^^^^^^(1) j _ ^^^^ . (^(Vo5(^) • Vo)Vo5(2) j ^ Q _ (28c/) 
(The proof is done by explicit verification). 

In special cases it is possible to construct the four potentials 5^^^) and {S^^^^)i 
by fixing the parameters Ai, A2, and fi2 suitably, thus fulfilling the constraints 
(27b) and (28d,e,f); in some cases the integration freedom (the curl of some vector 
in (27a), and the gradient of some scalar in (28a,b,c)) has to be used in order to 
fulfil the constraints, (see Buchert et al. 1993c). 
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4.3. Practical realization 

and comparison with numerical simulations 



In order to realize the presented solution for generic initial conditions the following 
procedure is appropriate. 

o — * 

Let us start with a Gaussian random density contrast field S{X) as initial condi- 
tion. It is most convenient to express this field in terms of its discrete Fourier 
transform (see, e.g., Bertschinger 1992). Henceforth, we denote its Fourier com- 
ponents by (5, the Fourier sums may extend over n wave vectors K = {Ki, K2, K3). 
From the density perturbation we construct the initial peculiar-velocity potential 
in space using FFT ("Fast Fourier Transform"), here written for the "flat" 
background: 

S{K) = -^^A-\K)6{K) . (29) 

From S{K) we are able to construct all displacement vectors which are generated 
by second derivatives of <S with respect to Lagrangian coordinates. This can be 
done by simply multiplying and summing the components Ki,K2, of all n wave 
vectors K in space. We first use the following property of Fourier transforms: 

VoS{X){K) ^iK S{K) ; (30a) 

{S,i,k){X){K) ^ -KiKk S{K) . (306) 

Transforming the gradient back to the X— space we obtain the displacement vector 
of the first order perturbation. Transforming all 6 different Fourier components of 
the symmetric tensor gradient (5 i,fc) back to the X— space we are able to construct 
the three principal scalar invariants of it: 

I{S,i,k) = «5,i,i + 5,2,2 + «5,3,3 , (31a) 

II{S,i,k) = <5,l,l»S,2,2 + «5,i,i<S,3,3 + <S,2,2»5,3,3 — <S,i,2<S,2,l — <S,i,3<S,3,i — <S,2,3<5,3,2 (316) 
III{S^i^k) = «5,3,35,2,2«5,i,i - 5,3,35,2,15,1,2 - 5 3,25,2,35,1,1 + 

'5,3,25,2,15 1,3 + 53,15,2,35 1,2 — 5,3,152,25 1,3 . (31c) 

From the invariants (31b) and (31c) we construct the potentials 5^^) and S^^"-^ by 
using again FFT in analogy to (29). Note, that for the displacement vectors in 
(24) we need the gradients of these potentials which can be constructed again by 
multiplying with iK in K— space. 
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To construct the interaction terms, we have to use the components of WqS^'^^ and 
mix them up with the components of VoS^^\ This is, of course, only possible after 
the construction of S^'^^: After building the mixed invariants, for the longitudinal 
interaction term: 

ttL (--(S) . <-.(2) (-.(2) (-.(2) 0(2) <-• 0(2) 

^^mix - '^,1,1'^,2,2 + '^,l,l'^,3,3 + '^,2,2'^,3,3 - '^,l,2'^,2,l " '^,1,3'>,3,1 " '^,2,3'^,3,2 

and for the transverse interaction terms: 

i,j,k= 1,2,3 cyclic ordering , {31e, f,g) 

we are able to construct the four generating functions of the interaction terms S^^^^ 
and 5*^"^'^^ via FFT, and from those functions the components of the displacement 
vectors. 

The presented model is currently being compared with numerical simulations. 
Two of the comparisons (Buchert et al. 1993a, b) concern the cross-correlation of 
generic density fields as predicted by the Lagrangian perturbation solutions with 
N-body simulations. This is done for various different power spectra similar to the 
work by Coles et al. (1993). Another comparison (Buchert et al. 1993c) performs 
a study of special cluster models using a tree-code. The purpose of the latter work 
is to learn about the principal effects of the different orders at high resolution. 
In this study the local forms discussed in 4.2 can be used, the models employed 
are examples of third-order solutions which are fully expressible in closed form 
without using a numerical Poisson solver. 

4.4. Discussion 

In the sequel all statements are understood in terms of the restriction (5). Similar 
statements will hold for the general case. 

The infinite sequence of a perturbation ansatz of the form: 

00 

/ = aX + 5^£y^) (32) 

£=1 
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contains in its spatial part integrals of invariants which involve sums and products 
of first derivatives of the perturbation fields pf. All these invariants are at most 
cubic in the products. This fact implies that the sequence (32) starts with those 
terms listed in (24) which grow as g^, q^z and q^^^, the whole rest of the sequence 
is made up of interaction terms resulting from products of perturbations which 
are not directly expressible in terms of the initial conditions. In view of this it is 
natural to retain those three terms as "main body" of the perturbation sequence 
for modeling the early non-linear regime, thus truncating all interaction terms. 
This way a truncated model is obtained which is the simplest third order model 
for the study of large-scale structure. In the sequel we list the advantages of the 
truncated model /*™"'^ = aX + q^VoS^^^ + q^^^oS^'^^ + q^^^VoS^^"^: 

• The model respects all three scalar invariants of the initial displacement 
tensor (5i,fc). 

• This (and only this) model maps the intial condition to the final stage 
directly without 

the need to construct the displacement vectors iteratively. 

• The model is purely longitudinal. 

However, the interaction terms (not contained in grow at the rates g^^^ 

and q^^g comparable with g"^^. The comparison with numerical simulations will 
show, whether neglection of these terms is meaningful. Note, that a similar growth 
rate does not imply similar importance for the clustering process, because that 
depends on the detailed spatial structure of the source terms in (24i-l) compared to 
the source term in (24h). The other question, whether the local forms discussed in 
section 4.2 might be useful also for generic initial conditions will also be resumed 
in future work. 

We finally emphasize that we can expect the Lagrangian perturbation ansatz to 
apply until a{t)a{to) = 0{1), where a{to) = §£ is the r.m.s. amplitude of the initial 
density contrast field. At this "accumulation point" all orders of the approxima- 
tions result in displacements of the same order as the first order displacements. 
Obviously, this restriction is much less severe as the bound \S\ < 1 in the Eulerian 
perturbation theory, since the density itself is not bounded in the Lagrangian 
framework. 
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Note: A REDUCE program of all equations in this paper and the solution can be 
obtained from the author via electronic mail: TOB at ibma.ipp-garching.mpg.de. 
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APPENDIX A 



In this APPENDIX the constraint equations resulting from the equations (2a,b,c) 
are Usted. Since we restrict ourselves to irrotational flows, a LEMMA proved in 
BE93 applies, which enables us to reduce the integrability conditions (2a,b,c) 
for the irrotationality of the gravitational fleld strength to the conditions (BE93: 
5d,e,f) for the irrotationality of the velocity. The components of the vorticity 
vector in Eulerian space uj = ^Vx x v, v = f, are expressed in Lagrangian space 
and are listed up to the second order for the ansatz (3) in (BE93, APPENDIX). 
The additional third order terms resulting from the same ansatz read explicitly 
(i,j,k = 1,2,3, cyclic ordering): 

The second order terms do not imply restrictions in the case of solutions of the 
form (6), since a functional relationship between the peculiar-acceleration and the 
peculiar- velocity is assumed (see BE93). The third order terms depend on the 
solutions (6a) and (6b). We insert them into the constraint equations and obtain: 

The terms which are of third order in the first order perturbations cancel for the 
present ansatz. We are left with constraints which result from interaction terms 
of first and second order perturbations. We conclude that interaction of the 
perturbations add trans versality to the third order solution which is not covered 
by the longitudinal ansatz. The Wronskian qzzqz — qzqzz is non-zero, because qz{t) 
and qzz{t) are linearly independent solutions. Thus, in order to fulfill the constraint 
equations for the longitudinal ansatz, we have to assure a functional relationship 
of the form Vo>S(2) = W(Vo>S(i)) with arbitrary W. Since AqS^^^ = 2/1(5^,^), 
we obtain the following equation for the functional W (a prime, here, denotes 
derivative with respect to the argument Vq^^^^; without loss of generality we can 
set 5(1) = Sto, see APPENDIX B): 

>V7(5 i,fc) - 2//(5 ,,fc)to = . 
Neglecting interaction terms altogether, the third order terms are not constrained. 
The constraints are removed by the transverse part derived in section 3.5.. 
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APPENDIX B 



In this APPENDIX a short discussion is given of the uniqueness of the solutions 
obtained. In particular, I shall demonstrate this in the case of the first order 
solution, and for the transverse part of the third order solution. 

For the realization of the first order solution we replace the perturbation potential 
S^^^ by the initial peculiar- velocity potential StQ. With this replacement the first 
order part of the solution (24) reduces to the "Zel'dovich approximation". This 
can be done without loss of generality as will be shown below: 
The solution to equation (24f) can be written as 

for any function ip which obeys the Laplace equation AqV' = 0. Thus, if we 
preserve the boundary conditions used for the initial condition S (here assumed 
to be periodic), then the only solution of the Laplace equation is V' = ^{t)- Without 
loss of generality we can set £{t) = 0. According to the form of the solution (24), 
this argument holds for any time. q.e.d. 

For the transverse part of the third order solution a vector potential S'^^'^^ has 
been introduced (eq. (19a)). On this a gauge condition has been imposed (eq. 
(19b)), which enables to write the transverse third order solution (eqs. (24j,k,l)) 
in terms of three Poisson equations. This is especially useful for technical reasons 
to realize the solution (4.3). In the following it is shown that the condition (19b) 
does not restrict the generality of the solution: 

Firstly, note that the vector S = — Vq x S'^^'^^ remains transverse, i.e., Vq ■ S = 0, if 
we add the gradient of an arbitrary scalar function O to the vector potential S'^^^K 
The divergence of this vector potential is yet unspecified. If we impose the gauge 
condition (19b), then we have to choose the gauge function fl in such a way that 

Aon = -Vo • .S^^'^) . 

This is always possible according to a theorem by Brelot on the existence of 
solutions of Poisson equations (see Friedman 1963). Imposing the gauge condition 
(19b) is sufficient to solve the equations (24j,k,l), since (19b) implies Vo(Vo-5(3'^)) = 
required for the derivation of these equations. q.e.d. 

A comprehensive discussion of the problem of uniqueness in general perturbation 
solutions will be given in a separate paper (Ehlers & Buchert 1993). 
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